Data Display

Default

plot(mtq)

Only geometry

plot(st_geometry(mtq))

Centroids Extraction

mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)

Distance Matrix

mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]

Units: [m]
          [,1]     [,2]      [,3]      [,4]      [,5]
[1,]     0.000 35297.56  3091.501 12131.617 17136.310
[2,] 35297.557     0.00 38332.602 25518.913 18605.249
[3,]  3091.501 38332.60     0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702     0.000  7177.011
[5,] 17136.310 18605.25 20226.198  7177.011     0.000

Polygons Aggregation

Simple union:

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")

Aggregation according to a grouping variable:

library(dplyr)
mtq_u2 <- mtq %>% 
  group_by(STATUT) %>% 
  summarize(P13_POP=sum(P13_POP))
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u2), add=T, lwd=2, border = "red", col=NA)

Buffer Zone

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")

Polygon Intersection

m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586), 
           c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)

Voronoi Polygons

google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)

mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join=st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col='lightblue')

Data Preparation

library(sf)
library(dplyr)
# Import geo layers
## Communes of Seine Maritime
sm <- st_read(dsn = "data/dep76/seine_maritime.geojson", stringsAsFactors = F, quiet=TRUE)
## French departements
dep <- st_read(dsn = "data/dep76/dep.geojson", stringsAsFactors = F, quiet=TRUE)

# change projection (lambert93)
sm <- st_transform(sm, 2154)
dep <- st_transform(dep, 2154)

# Import dataset  
csp <- read.csv("data/dep76/data_seine_maritime.csv")
# merge geolayer and dataset
sm <- merge(sm, csp, by="INSEE_COM", all.x=TRUE)

# Prepare some additional contextual informations
# Extract label of main cities
cities = c("Rouen","Fécamp","Le Havre","Dieppe","Le Tréport")
labels = sm %>% filter(LIBELLE %in% cities)

library(osmdata)
# Get major roads from osm
bb      <- sm %>% st_transform(4326) %>% st_bbox()
q       <- opq(bbox = bb,timeout = 180)
qm      <- add_osm_feature (q, key = 'highway',value = 'motorway',value_exact = FALSE)
qt      <- add_osm_feature (q, key = 'highway',value = 'trunk',value_exact = FALSE)
qp      <- add_osm_feature (q, key = 'highway',value = 'primary',value_exact = FALSE)

motorway<- osmdata_sf(qm)
trunk   <- osmdata_sf(qt)
primary <- osmdata_sf(qp)

roads    <- c(primary,trunk,motorway)$osm_lines %>% st_transform(st_crs(sm))
roads.geom = st_intersection(st_geometry(roads),sm) 

# Get the shape of the main river "La seine" 
qr    <- q %>% add_osm_feature (key = 'waterway') %>% add_osm_feature(key="name:fr",value="La Seine")
river <- osmdata_sf(qr)
river.geom <- st_geometry(river$osm_lines %>% filter(name.fr=="La Seine")) %>%
  st_transform(st_crs(sm))

Proportional symbols

# (Very) simple map
library(cartography)
plot(st_geometry(sm))
propSymbolsLayer(sm, var = "act")
title("Active Population")

# Custom map of active population
par(mar=c(0.2,0.2,1.4,0.2))
bb <- st_bbox(sm)
# the bbox is used to center the map on the Seine Maritime depatement
plot(st_geometry(dep), col = "ivory", border="ivory3",  bg="azure", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
plot(st_geometry(sm), col="cornsilk2", border = NA, lwd = 0.5, add=T)
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
plot(st_geometry(river.geom),col="azure",lwd = 3,add=TRUE)
propSymbolsLayer(sm, var = "act", col="darkblue", inches = 0.6, 
                 border = "white", lwd=0.7, symbols = "square",
                 legend.style = "e", legend.pos="topleft",
                 legend.title.txt = "Labor Force\n(2014)", 
                 legend.values.rnd = 0)
labelLayer(labels,txt="LIBELLE",halo=TRUE,col="white",bg="black",cex=1) 
# Scale Bar
barscale(size = 10)
# North Arrow
north(pos = "topright", col = "darkblue")
# Full layout
layoutLayer(title = "Workforce in Seine-Maritime", 
            sources = "Insee, 2018, OpenStreetMap contributors",
            author = "Kim, Tim & Comeetie, 2018", 
            col = "darkblue", coltitle = "white", tabtitle = TRUE, 
            frame = TRUE, scale = NULL, north = FALSE)

Typology Map

# To display qualitative data
# modalities
mod <- c("agr", "art", "cad", "int", "emp", "ouv")
# labels in the legedn
modlab <- c("Agriculteurs", "Artisans","Cadres", "Prof. Inter.", "Employés", "Ouvriers")
# colors
cols <- c("#e3b4a2", "#a2d5d6", "#debbd4", "#b5dab6", "#afc2e3", "#e9e2c1")

par(mar=c(0.2,0.2,1.4,0.2))
plot(st_geometry(dep), col = "ivory", border="ivory3", bg="azure",
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
typoLayer(sm, var = "cat", 
          border = "ivory", lwd = 0.5, 
          legend.values.order = mod,
          col = cols,
          add=TRUE, legend.pos = "n")
plot(st_geometry(river.geom),col="azure",lwd = 3,add=TRUE)
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
# functions are dedicated to legend display
labelLayer(labels,txt="LIBELLE",halo=TRUE,col="white",bg="black",cex=1) 
legendTypo(title.txt = "Dominant Socio-Professional\nCategory", 
           col = cols, 
           categ = modlab, 
           nodata = F)
barscale(size = 10)
north(pos = "topright", col = "darkblue")
layoutLayer(title = "Workforce Distribution in Seine-Maritime", 
            sources = "Insee, 2018, OpenStreetMap contributors",
            author = "Kim, Tim & Comeetie, 2018", 
            col = "darkblue", coltitle = "white", tabtitle = TRUE,
            frame = TRUE, scale = NULL, north = FALSE)

Choropleth Map

# Compute the share of "managers" in the active population
sm$pcad <- 100 * sm$cad / sm$act

# The getBreaks() function is used to classify the variable
bks <- getBreaks(v = sm$pcad, method = "quantile", nclass = 6)

# The carto.pal() function give access to various cartographic color palettes
cols <- carto.pal("green.pal", 3,"wine.pal",3)

# Create the map
par(mar=c(0.2,0.2,1.4,0.2))
plot(st_geometry(dep), col = "ivory", border="ivory3", bg="azure",
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
choroLayer(sm, var = "pcad", breaks = bks, 
           col = cols, border = "grey80", 
           legend.values.rnd = 1,
           lwd = 0.4, legend.pos = "topleft", 
           legend.title.txt = "Share of managers (%)", add=T)

plot(st_geometry(river.geom),col="azure",lwd = 3,add=TRUE)
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
# functions are dedicated to legend display
labelLayer(labels,txt="LIBELLE",halo=TRUE,col="white",bg="black",cex=1) 
# Add a layout
layoutLayer(title = "Managers", 
            sources = "Insee, 2018, OpenStreetMap contributors",
            author = "Kim, Tim & Comeetie, 2018", 
            theme = "green.pal", 
            col = "darkred", coltitle = "white", 
            tabtitle = TRUE, 
            frame = TRUE, scale = 10)
north(pos = "topright")

We could create the same map on a cartogram based on the active population stock.

library(cartogram)
sm_c1 <- cartogram_cont(sm, "act", 3)
par(mar=c(0.2,0.2,1.4,0.2))
plot(st_geometry(dep), col = "ivory", border="ivory3", bg="azure", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
choroLayer(sm_c1, var = "pcad", breaks = bks, 
           col = cols, border = "grey80",            
           legend.values.rnd = 1,
           lwd = 0.4, legend.pos = "topleft", 
           legend.title.txt = "Share of managers (%)", add=T)
# functions are dedicated to legend display
labelLayer(sm_c1 %>% filter(LIBELLE %in% cities), txt="LIBELLE", halo=TRUE, col="white",
           bg="black", cex=1) 
# Add a layout
layoutLayer(title = "Managers", 
            sources = "Insee, 2018",
            author = "Kim, Tim & Comeetie, 2018", 
            theme = "green.pal", 
            col = "darkred", coltitle = "white", 
            tabtitle = TRUE, 
            frame = TRUE, scale = 10)
north(pos = "topright", south = TRUE)

Gridded Map

These maps may allow to hide or, at least, diminish the MAUP.

# Create a grid based on sm (cell area = 4 square km)
grid <- getGridLayer(x = sm, cellsize = 4000*4000, 
                     type = "regular", var = c('cad', 'act'))
# Compute the share of managers
grid$pcad <- 100 * grid$cad / grid$act
# Display the map as choropleth layer
par(mar=c(0.2,0.2,1.4,0.2))
plot(st_geometry(dep), col = "ivory", border="ivory3", bg="azure", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
choroLayer(grid, var = "pcad", breaks=bks, 
           col = cols, border = "grey80", 
           lwd = 0.4, legend.pos = "topleft", 
           legend.title.txt = "Share of managers\n(en %)", add=T)
plot(st_geometry(river.geom),col="azure",lwd = 3,add=TRUE)
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
# functions are dedicated to legend display
labelLayer(labels,txt="LIBELLE",halo=TRUE,col="white",bg="black",cex=1) 
layoutLayer(title = "Managers", 
            sources = "Insee, 2018, OpenStreetMap contributors",
            author = "Kim, Tim & Comeetie, 2018", 
            theme = "green.pal", 
            col = "darkred", coltitle = "white", 
            tabtitle = TRUE,
            frame = TRUE, scale = 10)
north(pos = "topright")

It is also possible to create hexagonal grids.

# Create a grid based on sm (cell area = 4 square km)
grid2 <- getGridLayer(x = sm, cellsize = 4000*4000, 
                     type = "hexagonal", var = c('cad', 'act'))
# Compute the share of managers
grid2$pcad <- 100 * grid2$cad / grid2$act
# Display the map as choropleth layer
par(mar=c(0.2,0.2,1.4,0.2))
plot(st_geometry(dep), col = "ivory", border="ivory3", bg="azure", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
choroLayer(grid2, var = "pcad", breaks=bks, 
           col = cols, border = "grey80", 
           lwd = 0.4, legend.pos = "topleft", 
           legend.title.txt = "Share of managers\n(en %)", add=T)
plot(st_geometry(river.geom),col="azure",lwd = 3,add=TRUE)
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
# functions are dedicated to legend display
labelLayer(labels,txt="LIBELLE",halo=TRUE,col="white",bg="black",cex=1) 
layoutLayer(title = "Managers", 
            sources = "Insee, 2018, OpenStreetMap contributors",
            author = "Kim, Tim & Comeetie, 2018", 
            theme = "green.pal", 
            col = "darkred", coltitle = "white", 
            tabtitle = TRUE,
            frame = TRUE, scale = 10)
north(pos = "topright")

Smoothed Map

smoothLayer() uses functions from package SpatialPosition to compute Stewart’s potentials of population.
The computation of potentials could be considered as a spatial interpolation method such as inverse distance weighted interpolation (IDW) or kernel density estimator. These models aim to estimate unknown values of non-observed points from known values given by measured points. Cartographically speaking, they are often used to get a continuous surface from a set of discrete points. However, Stewart model is mainly a spatial interaction modeling approach, with a possible secondary use for spatial interpolation.

grid$cad100 <- grid$cad * 100
par(mar=c(0.2,0.2,1.4,0.2), bg="azure")
plot(st_geometry(dep), col = "ivory", border="ivory3", 
     xlim = bb[c(1,3)], ylim =  bb[c(2,4)])
smoothLayer(x = grid, var = "cad100", var2 = "act", 
            typefct = "exponential", 
            span = 4000, beta = 2, breaks = bks, col = cols,
            legend.pos = "topleft", mask = st_buffer(sm, 0),
            legend.values.rnd = 1,
            legend.title.txt = "Share of managers* (%)", 
            border = "grey90", lwd = 0.2, add=T)
plot(st_geometry(river.geom),col="azure",lwd = 3,add=TRUE)
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
# functions are dedicated to legend display
labelLayer(labels,txt="LIBELLE",halo=TRUE,col="white",bg="black",cex=1) 
layoutLayer(title = "Managers", 
            sources = "Insee, 2018, OpenStreetMap contributors",
            author = "Kim, Tim & Comeetie, 2018", 
            theme = "green.pal", 
            col = "darkred", coltitle = "white", 
            postitle = "center",
            frame = TRUE, scale = 10)
north(pos = "topright")
text(x = 488000, y = 6921000,adj = 0,
     font = 3, cex = 0.8,
     labels = "* Potential smoothing\n  exponential function\n  span = 4 km, beta = 2")

Cheat Sheet

The cheat sheet displays a quick overview of cartography’s main features:

ggplot2

Proportional symbols layer

library(ggplot2)

labgg = labels %>% st_centroid() 
labgg$x = st_coordinates(labgg)[,1]
labgg$y = st_coordinates(labgg)[,2]

ggplot() +
  geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
  geom_sf(data = river.geom, colour = "azure",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  geom_sf(data = sm %>%  st_centroid(),
          aes(size= act), colour="#E84923CC", show.legend = 'point') +
  geom_text(data=labgg,aes(x=x,y=y,label=LIBELLE))+
  scale_size(name = "Active population",
             breaks = c(100,1000,10000),
             range = c(0,20)) +
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(sm)[c(1,3)],
           ylim = st_bbox(sm)[c(2,4)]
  ) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(title = "Active population",
       caption = "Insee, 2018\nKim & Tim, 2018",x="",y="")

Choropleth layer

ggplot() +
  geom_sf(data = dep, colour = "ivory3",fill = "ivory") +
  geom_sf(data = sm, aes(fill = pcad), colour = "grey80") +
  geom_sf(data = river.geom, colour = "azure",size=2) +
  geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
  geom_text(data=labgg,aes(x=x,y=y,label=LIBELLE))+
  scale_fill_gradientn(name = "Share of managers (%)",
            
                       colours = carto.pal("green.pal", 3,"wine.pal",3),
                       values=bks/max(bks))+
  coord_sf(crs = 2154, datum = NA,
           xlim = st_bbox(sm)[c(1,3)],
           ylim = st_bbox(sm)[c(2,4)]
  ) +
  theme_minimal() +
  theme(panel.background = element_rect(fill = "azure",color=NA)) +
  labs(title = "Managers",
       caption = "Insee, 2018\nKim & Tim, 2018",x="",y="")

tmap

TODO Timothée

Proportional symbols layer

These interactive maps, as appealing as they seem to be, are not really suitable for presenting geostatistical information. For porportional symbols both leaflet and mapview lack a proper way to build lengends. Nonetheless, they can be really useful for exploratory data analysis. Both librairies are quite similar with some advantage for mapview (e.g the possibility to render to canvas for bigger datasets) and some other for leaflet (e.g using a customized projection different from web-meractor)

library(mapview)
# convert polygon to point and compute the circle radius for porportional area
actifs = sm %>% st_centroid() %>% mutate(radius = sqrt(act)) %>% arrange(desc(act))
# build a map by providing a map tile provider and some options for the circles.
# with mapview the circle size will stay constant whatever the zoom level is.
mapview(actifs, map.types = "Stamen.TonerLite", cex="radius", legend=FALSE,
        col.regions="#217844", lwd=0, alpha=0.4)

With leaflet data must be provided in long/lat (st_transform(4326)) and will be converted by default to the webmercator coordinate reference system (CRS). You may choose between circles with radius specified in meters (addCircle) or in pixels (addCircleMarkers). If specified in pixel, the radius will stay the same whatever the zoom level is. If specified in meters, it will evolve with zoom level.

library(leaflet)
leaflet(actifs %>% st_transform(4326)) %>% addProviderTiles(providers$Stamen.TonerLite) %>%
  addCircles(radius = ~radius * 30, popup = ~paste0(LIBELLE," : ",act,"actifs"),
             fillColor = "#217844", stroke = FALSE, fillOpacity = 0.4)

Gridded Map

mapview and leaflet enable the use and styling of sf polygons. We may for example reuse the grid built previously to show the share of managers in Seine-Maritime. Default choropleth can be obtained with just three arguments with mapview.

mapview(grid,map.types = "Stamen.TonerLite",zcol="pcad")

leaflet enables fine tuning of the color scale and legend.

# define a color scale function from colors and breaks
color.scale <- colorBin(cols, domain = range(bks), bins=bks)

leaflet(grid  %>% st_transform(4326)) %>% addProviderTiles(providers$Stamen.TonerLite) %>%
  addPolygons(fillColor = ~color.scale(pcad), stroke = FALSE, fillOpacity = 0.7) %>% 
  addLegend(colors = cols, labels = round(bks[1:(length(bks)-1)]*10)/10,
            title = "Share of Manager (%)")

Projection

leaflet enables to use other projections than the web-mercator projection for webmaps. You may deal with tiles built with another projection system such as Lambert 93 by defining a custom leaflet crs (see gis.stackexchange resolution-from-wmts-getcapacilities-scaledenomin and mathematics-behind-converting-scale-to-resolution to deal with WMTS tiles). We use this approach here to deal with Hypsometric tints tiles (tiles colors which encode elevation) and the contour lines of chamois and ibex population area provided by the french national office of hunting and wild animals.

epsg2154 <- leafletCRS(crsClass = "L.Proj.CRS", code = "EPSG:2154",
                       proj4def = "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3
                       +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
                       +units=m +no_defs",
                       resolution = 4838095.23807/cumprod(c(1,rep(2,15)))*0.00028,
                       origin= c(-357823.2365,7230727.3772)
)

tile_ortho <- "http://tiles.craig.fr/ortho/service?service=WMTS&request=GetTile&version=1.0.0&
layer=ortho_2016&style=normal&tilematrixset=lambert93&format=image%2Fjpeg&tilematrix={z}
&tilerow={y}&tilecol={x}"

tile_mnt <- "http://tiles.craig.fr/mnt/service?service=WMTS&request=GetTile&version=1.0.0&
layer=relief&style=default&tilematrixset=lambert93&format=image%2Fpng&tilematrix={z}&
tilerow={y}&tilecol={x}"

tile_attrib <- "Map data &copy;
<a href='https://www.craig.fr/contenu/1377-flux-tuile-wmswmtstms'> craig / IGN </a>
and <a href='http://www.oncfs.gouv.fr/Cartographie-ru4/Le-portail-
cartographique-de-donnees-ar291'>ONCF</a>"

bouq = read_sf("./data/BOQ_2017_massif_uni/BOQ_2017_massif_uni_L93.shp")
cha = read_sf("./data/CHA_2017_massif_uni/CHA_2017_massif_uni_L93.shp")

leaflet(options = leafletOptions(worldCopyJump = F, crs = epsg2154)) %>%
 # addTiles(urlTemplate = tile_ortho,attribution = tile_attrib) %>% 
  addTiles(urlTemplate = tile_mnt,attribution = tile_attrib) %>% 
  addPolygons(data = bouq %>% st_transform(4326), fillColor = '#50162D',color = '#772A7F',
              fillOpacity = 0.05,opacity = 1, weight = 2.5, dashArray = "3") %>% 
  addPolygons(data = cha %>% st_transform(4326), fillColor = '#88AA00',color = '#CCFF00',
              fillOpacity = 0.05,opacity = 1, weight=2.5, dashArray = "3") %>% 
  setView(2.692174, 45.067230, 4)